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Abstract 

This paper examines an inverse problem in thermal imaging, that of recovering a void in a 
material from its surface temperature response to external heating. Uniqueness and contin- 
uous dependence results for the inverse problem are demonstrated and a numerical method 
for its solution developed. This method is based on an optimization approach, coupled with 
a boundary integral equation formulation of the forward heat conduction problem. Some 
convergence results for the method are proved and several examples are presented using 
computationally generated data. 
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1 Introduction 


Thermal imaging is a technique of recent interest for the nondestructive evaluation of mate- 
rials. This method attempts to characterize the internal structure of a sample (perhaps to 
locate flaws-disbonds, bubbles, corrosion, etc.) by using its surface temperature response to 
an external heating. Some recent work and techniques in this subject are detailed in [3], [4], 
[5], [7] and [9]. 

In this paper the problem of detecting and identifying an unknown internal void in a 
planar domain using thermal imaging is examined. The void could represent a defect in the 
material, or it could be a feature which is supposed to be present, e.g., a conduit, whose 
location or geometry is to be assessed. The focus is on the case in which the thermal 
stimulus, an applied heat flux at the boundary of the sample, is periodic. After separating 
the temporal and spatial variables, one obtains an inverse or domain identification problem 
for an elliptic equation. Results concerning the uniqueness and continuous dependence for 
the inverse problem will be examined and an algorithm for the numerical recovery of the 
void will be presented. This algorithm will be applied to examples using computationally 
generated data with and without noise. 

The outline of the paper is as follows. Section 2 concerns the mathematical formulation 
of the forward heat conduction problem with periodic heating and demonstrates how this 
leads to an inverse problem for an elliptic equation. Section 3 contains an identification 
result for the inverse problem and shows that an internal void is determined by the bound- 
ary temperature response to external heating. This section also contains results concerning 
sensitivity or continuous dependence for estimates of the internal void based on the bound- 
ary measurements. These results rely on reformulating the heat conduction problem as an 
integral equation on the boundary of the sample. Section 4 examines a numerical method 
for the solution of the forward problem based on the boundary integral formulation and its 
incorporation into a least-squares routine for solution of the inverse problem. It is shown 
that with reasonable hypotheses on the class of voids, the numerical method will converge 
to the solution of the least-squares formulation of the inverse problem. Section 5 presents 
the results of this algorithm applied to computationally generated data. 
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2 Mathematical Formulation 


The sample (without void) to be tested will be denoted by f l, a bounded region in IR 2 with 
C 2 boundary. The internal void will be denoted by D , where D CC with C 2 boundary. 
The function T(t, x ) will denote the solution to the heat equation 

0 in ft \ D 

g(x,t) 

0 

T 0 (x) 

where v is the outward unit normal vector field on the boundary of fl \ D, k is the thermal 
diffusivity of f), a is the thermal conductivity of f 1, T 0 denotes the initial temperature of 
the region and g is the heat flux at the boundary. Both n and a are assumed to be known 
constants. Of course it is assumed that g is not identically zero. 

The heat flux g(x,t ) will be assumed to be periodic in time with known frequency ^ so 
that 

g(x,t) = Re{e ,u "#(*)} 

for some complex valued function g(x). Actually, g would also generally include a steady- 
state term as well, but since only the periodic response is of interest, this term can be ignored. 
Under this assumption one can separate variables to find that T(x,t) = Re{T(x)e twi } where 
T(x) satisfies 

= 0 in ft \ D 

= g(x) (2.1) 

= o, 

at least for time large enough so that the initial conditions do not matter. Note that the 
function T(x) solving (2.1) will be complex- valued, consisting of a real, or in phase, and 
imaginary, or out of phase part. 
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The inverse problem of interest is the following: Given a known applied heat flux < 7 , can 
the shape and location of the void D be uniquely determined from measurements of T on 
the boundary of 0? Provided a uniqueness result holds, one would also like to know whether 
D depends continuously on measurements of T on <9fi, that is, how sensitive estimates of 
D are to noise in the data. Finally, one would like an efficient computational algorithm for 
recovering an estimate of D from actual data. 


3 The Inverse Problem 

3.1 Uniqueness 

A uniqueness result for the inverse problem follows easily from basic facts about elliptic 
operators. 

Theorem 3.1 ( Uniqueness ) Let D\ and D 2 be two subdomains of and T\ and T 2 the 
corresponding solutions to equation (2.1) with nonzero Neumann data g. Let S be a portion 
of dLt with positive measure. Then T\ = T 2 on S implies that D\ = D 2 . 

Proof: The functions T\ and T 2 have the same Neumann data g and, since Ti and T 2 
agree on S , the same Cauchy data. Unique continuation for elliptic operators implies that 
Tr and T 2 agree on fi\ {D x U D 2 ). Then, for example, the function T 2 has a vanishing normal 
derivative on the region D\ \ D 2 , hence is constant on D\ \ D 2 . If D\ \ D 2 ^ 0 then by the 
maximum principle T 2 must be constant throughout f 1 \ D 2 , contradicting g ^ 0. A similar 
argument shows that if D 2 \ D\ ^ 0 then T\ would be constant on Ll\D\, again contradicting 
g 0, hence D\ = D 2 . 

3.2 Boundary Integral Formulation 

In order to investigate continuous dependence and numerical methods for the recovery of 
D it will convenient to reformulate the heat conduction problem as a boundary integral 
equation. This offers the advantage that one only has to solve for the temperature on the 
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boundary of the sample, rather than the interior. Since the boundary is the only place the 
temperature is measured, this is the only place its value is needed. 

Use r(r) to denote a fundamental solution or Green’s function for the operator (A - ^). 
Such a fundamental solution is given by 

f(r) = 

= (ker( r \/w/K) + iVe\{ryJu/K)) 

where K 0 is the zero order modified Bessel function of the second kind and ker and kei denote 
the kelvin functions. Efficient routines for computing the kelvin functions can be found in 
[1]. Define 

r(*,y) = f(|ar-y|). 

The function T satisfies the heat equation in the y variable for fixed x, except at x = y, 
where it has a logarithmic singularity. Standard potential theory arguments (see [6], chapter 
3) show that the elliptic problem given by equation (2.1) can be formulated as a boundary 
integral equation, 

“ \ T{x) + 4)\ D) T{y)d "> T(x ' v)iS ’ = 4>y» r(x ’ *Wv)iS, < 31 ) 

for each x 6 <9(0 \ D ) where d Uy is the normal derivative in the y variable and dS y is surface 
measure. We will use K(x,y) to denote the kernel d Uy Y{x,y) for x,y G d(Yl \ D) and use S 
to denote the operator 

{S<t>)(x) = f K (x, y)<t>{y) dS y . (3.2) 

The operator S is bounded and compact on C(d(fl \ D )), the space of continuous functions 
on \ Z)), hence — ~I + S is a second kind Fredholm operator. Uniqueness of the solution 
to equation (3.1) follows from uniqueness of the solution to the forward problem. By the 
Fredholm alternative, the equation ( — + S)<f) = g is solvable, at least for smooth enough 
g. In particular, ( — ^ + 5*) 1 exists and is bounded on C(d(fl\D)). The solution to equation 
(3.1) yields the temperature T(x ) on dfl and dD. If needed, the temperature for x £ Q \ D 
can be found from the relation (Green’s third identity) 

T ( x ) = - r ( x ,y) 9 (y))dS y . 
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3.3 Restrictions on the Domain 

In order to obtain continuous dependence results, a few restrictions on the class of voids 
and their parameterization are needed. First, let us use C 2 [ 0, 1] to denote the space of C 2 
functions on [0, 1] where the endpoints 0 and 1 are identified with each other. This space 
can be normed by \\<f>\\ = sup t€[0>1] \D p <j)\, \0\ < 2. We assume that D depends on finitely 
many parameters, D = D(q) with q E Q CC IR m and where: 

(a) qi = q 2 implies D(qi) = D(q 2 ) (unique parameterization). 

(b) D(q) C f Y CC for q E Q ( D(q ) stays away from dtt). 

(c) The closed curves dD(q) are parameterized as x(q,t) = (xi(q,t),x 2 (q,t)) for q € Q, 
0 < t < 1, with X{(q,t ) a C 2 function of t for each q E Q and ^ = \f{x\) 2 + {x^) 2 
bounded away from zero. Also, the map q —* x{q , t) is continuous from IR to C 2 \ 0, 1]. 
Ba^ed on the above assumptions it is not difficult to show that: 

(d) For each x £ dD(q ) there exists an open ball, B c (x ), of radius e around x with 
t independent of x and q such that the curve dD(q) (1 B c (x) is parameterized by 
(xi(q,t), x 2 (q, t)) with t in some connected interval. 

The continuity of q -> C 2 and compactness of Q imply that the C 2 norm of x(q,t) as a 
function of t is uniformly bounded over q € Q- Also, for q € Q the families of functions 
x(q,t ), x'(q,t ) and x"(q,t) are equicontinuous in t,, where the prime denotes differentiation 
with respect to t. 

Lemma 3.1 The family of functions {K(x(q,s),x(q,t))- q € Q}, are (uniformly) equicon- 
tinuous in s and t , that is, for each t > 0 there is a 5 > 0 so that 

\K(x(q,s),x(q,t))- K(x{q,s'),x(q,t'))\ < e 

for all ( s,t ) and (s'jt 1 ) with 

\s - s’ I <6, 1 1- a I < 8 , 

and 8 does not depend on s, t or q. 
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Proof: For brevity let us suppress the dependence of z() and A"() on q and write simply 
A' CM) for I<(x(s),x(t)). Also, since T(x,y) = log(|x - y\) + G{x,y ) where G is smooth in 
x and y, we will prove the lemma assuming that r(z,j/) = log(|z — j/|); it is easy to check 
that the smooth term makes no difference in the proof. 

The stated regularity for K holds on the compact set {( s,t ) £ [0, 1] x [0, 1], |s — 1| > e, 
q € Q}, where e is any number greater than zero, for on this set K is at least C 2 . We need 
only to show that K ($, t ) is uniformly continuous on the set (s, /) £ [0, 1] x [0, 1], |.s — 1 1 < e. 

Suppose the boundary of D near a point Xo is parameterized by x(t) = x 2 (t)). By 

making an appropriate translation and rotation it may be assumed that x 0 = z(0) = (0,0) 
and that the boundary is oriented so that the unit normal outward vector in (z 5 , £ 2 ) coor- 
dinates is (0, 1). In this case Taylor’s theorem can be used to expand x{t) as 

x i(t) = a o t + — t 2 + R](t) 

X 2(t) — yt 2 ^2 (t) 

where 


a 0 = x\ (0) 

«i = x"(0) 

b = ^(O) 

«,(<) = i«(c) - x [ l( 0 ))t\ i = 1,2, 

and c is some point between 0 and t. The functions R 1 and R 2 are functions whose G 2 norms 
can be bounded in terms of the norms of x^ and x 2 . The unit normal vectors satisfy 

= * a (0 = bt + A ' 2 (0 
dS 

1/2 ^~dt ~ -x i(0 = -«° - ai ^ - ^(t). 

The kernel K(s,t ) is then given by 


dS_ _ 1 (xi(s) - xi + (x 2 (s) - x 2 (t))i/ 2 (0) dS 

dt 2n (zi(s) — xi(t)) 2 + (^ 2 ( 5 ) — x 2 {t')) 2 dt 

_ 1 (zips) - x{{t))^ 2 (t) - (z 2 (-s) - ^(QjfUQ) 

2tt (zi(s) - x x {t)) 2 + (z 2 (s) - x 2 {t)) 2 
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Substituting the above expressions for xi,X 2 ,x\ and x' 2 gives for the numerator, after some 
simplification, 

rium = -^(s - t ) 2 + R 2 (t)(s - t) a 0 + y (5 + t) + [R^s) - R x (t)](bt + R! 2 (t)) 

— (Ri( s ) ~ R‘i (t)) [«o + (i\t + R[(t)] — — (s 2 ~ 1 ( 0 • 

The denominator becomes 2i r times 

■ -| 2 '7 "2 

a 0 (t “ 5) + — (S 2 — t 2 ) + — /?i(£) + “(s 2 — ^ 2 ) + /?2(0 

Taylor’s theorem implies that 

R,(s) = R,(t) + (5 - <)#>(<) + i(> - 

for some point c between s and t. Substituting this into the expression for the numerator 
gives 

(s - t) 2 

num — - [— a^b + ( bt -f R2(t))Ri(c) — (&q + c i\t + R!^{ty) R f 2 {c) + a\R f 2 {t ) — bR\(t)] 

where c and c are between s and t . Doing the same for the denominator shows 

den ~ 2i r(s — t ) 2 ^{ciq T i?|(/)) 2 T (/? 2 (0) 2 T C\.s T 

where Ci and C 2 are functions of 5 and t which can be bounded in terms of the C 2 norm of 
x{t). Cancelling the common ( s — t ) 2 factor, the kernel can be written as 

, . _ 1 — apb + [bt + R f 2 (t)]R r ((c) — [aq — a\t — Ri(t)]R%(c) + a\R! 2 {t) — bR\(t ) 

' “ 4t r [ao + + [/^(f)] 2 + (7,5 + C 2 t 1 > 

Since /2j(0) = 0, the kernel is bounded through s = t, for the denominator is bounded 
away from zero for s and t in a sufficiently small neighborhood of zero. In fact, by Taylor’s 
theorem R\(t) = i?j(0) + R"(c)t = R"(c)t for some c between zero and t. In addition, since 
the functions C\ and C 2 can be bounded in terms of the C 2 norm of x(t ) (which is itself 
bounded uniformly over Q ), the denominator may be bounded uniformly away from zero 
for s and t in a neighborhood of 0, with the neighborhood and bound independent of xq 
and q. The families of functions x'j(t) and x"(t), as well as Rj(t) and R"(t ), are uniformly 
equicontinuous for q E Q, and hence so are both the numerator and denominator of equation 
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(3.4). Since the denominator of I\(s,t ) is bounded away from zero, it follows that the family 
of functions {I<(x(q, s), x(q, t))] q Q] is uniformly equicontinuous in s and t. ■ 

By parameterizing dD(q ) with x(t), 0 < t < 1 as above and parameterizing dD with x(t), 
\ <t <2, one can identify the boundary of d(D\D) with the fixed interval [0,2). A solution 
T(x ) to equation (3.1) can be identified with a function T(t ) on [0,2) by T(t ) = T(x(t)). 
For a function <f> defined on [0,2) let <^>i and (j) 2 denote the restriction of <f> to [0,1) and 
[1,2), respectively. We will work in the space of functions (f) for which <f > i is continuous and 
extends continuously to [0, 1] with <^i(0) = <^i(l) and for which (f> 2 is continuous and extends 
continuously to [1,2] with <f> 2 { 1) = <fo(2). We will denote this space by C[0,2] and norm it 
with the supremum norm. The solutions T{t) to (3.1) lie in this space. One can also identify 
the operator S (a function of q) with an integral operator on (7[0, 2]. Let us write I\(q,s,t ) 
instead of K(x(q,s),x(q,t)) and define 

2 J Q 

S(q)m = J o K(q,s,t)cj>{t)-^dt. (3.5) 

The same argument given in the proof of Lemma 3.1 shows that K(q,s,t) is uniformly 
equicontinuous for s and t in [1,2), i.e., for a;(s) and x(i ) on the boundary of Cl and q Q, 
(actually, K here is independent of q). For s G [0,1) and < € [1,2) (or vice-versa) K(q,s,t ) 
is C 2 in 5 and t , since in this case x(s) € dD , x{t ) € dfl and by assumption the boundaries 
do not intersect, so again K(q,s,t) is uniformly equicontinuous over q € Q- The kernel 
K will have a jump discontinuity at s = 1 or t = 1, since there x(t) jumps from dD to 
dfi. In summary, the family of functions {K(q, s, t ); q € Q}, is uniformly equicontinuous on 
[p,p+ 1) x [<?, q + 1) Pi q = 1)2, with simple jump discontinuities at s = 1 or t = 1. 

Another fact worth noting is that the map q — ► K ( q , s, f ) is continuous as a map from 
lR m to C[0, 2] x C[ 0, 2], that is, for any t > 0 there is a ^ so that \q - q\ <5 implies 

sup \I<(q,s,t) - K(q,s,t\ < e. 

*,<e[o,2] 

This follows directly from equation (3.4) and the fact that q — > x'j(q,t), q — * x'j(q,t), 
q _ > R'jit) and q — » R"(t) are all continuous as maps from IR m to the space of continu- 
ous functions. 
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3.4 Continuous Dependence 

Based on the above assumptions it is possible to prove a version of continuous dependence. 

Theorem 3.2 Let q n be a sequence in Q and T(q n ) the corresponding solution to equation 
(2.1) with D = D(q n ). Suppose T(q n ) — » T(q*) in C(d(Lt\ D )) for some q* E Q. Then 
q n -» <7*- 

Proof: The first task is to show that the map q — > S(q) is continuous. Let Si — S(q) denote 
the boundary integral operator for D = D{q) (considered as an operator on C[0,2j) and 
S -2 = S(q + Sq), where q E Q and Sq is some small perturbation in q. As remarked above, 
q — » K(q,s,t ) is continuous. It follows that q — >■ S(q) is also continuous as a map from lR’ n 
to the space of operators on C[0,2], for if \K(q + Sq,s,t) — K(q,s,t)\ < e for all 6 [0,2] 
then 

||(S , 2</>)(5) - (S'i<^)(s)|| < f \K(q + 8q,s,t) - I<(q,s,t)\\(f>(t)\dt 

JO 

< Ce \\ii 

for some constant C, that is, + 8q) — 5(^)|| < Ce for Sq sufficiently small. 

The next step is to show that q — » (/ — S(q))~ l is continuous. The operator (— |/ + .Si) 
is invertible and {—\l + Sf) can be inverted with a Neumann series as follows. First note 
that 

(-j/ + sy-’ = (-i/ + s, - (.?, - s 2 ))-' 

= (-\l + s,r'(r + (~\l + S,)(S, - Si))-\ 

Let R = — ( — j/ + S\)(Si — Si))- Given any e > 0 one can choose a |£g[ sufficiently small so 
that ||f?|| < e and, for e < 1, 

J-I + Sl )~\I + R)- 1 

"I + 5i) _1 (7 + R + R 2 ) 

~\l + S i)~' = R + R 2 + --- 


(~\l + S2)-' = ( 
= ( 

so that 
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and 


< 


( 3 . 6 ) 

( 3 . 7 ) 


K-jZ + s .)- 1 


(-^"i 


< 


m 

1 - \\ R \\ 

e 

1 - e’ 


Thus the map q — > ( — |/ + S(q))~ l is continuous and so the map q — > T(q) = ( — | + S(q))~ l g 
is continuous from IR m to ( 7 [ 0 , 2 ]. 

To complete the proof, suppose that the sequence q n does not converge to q *. Since <2 
is compact, some subsequence of qr n converges to some <7 € Q. It can be assumed that this 
is simply the sequence q n . However, continuity of q —> T(q ) means that T(q n ) — > T(q ), 
implying T(q*) = T(q) and contradicting the uniqueness Theorem 3.1. 


4 Numerical Methods 

4.1 Nystrom’s Method 

The following computational approach to the inverse problem is based on a least-squares 
formulation, finding the model void parameters which best fit the measured data by means 
of an optimization method. One drawback to this approach is the need to repeatedly solve 
the heat conduction problem (2.1). It is thus advantageous to have an efficient method for 
solving this equation. The boundary integral equation approach is such a method, and in 
this section we examine a technique for its solution. 

The boundary integral formulation ( 3 . 1 ) for the solution T to equation (2.1) can be 
written 

-\ T ( S ) + Jo K(s,t)T{t)dt = g(s) ( 4 . 1 ) 

where T(s) means T(x(s)) for the parameterization x(s ) of < 9 (fl \ D ). Surface measure 
has been included in the kernel K. Let tj and u>j j = 1 , ■ ■ • ,n denote the nodes and weights 
of a quadrature rule convergent on the space C[0,2], so that 

j=i 
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if /(<) is smooth enough. 
n, with 


Actually, we will consider a family of quadrature rules, indexed by 

lim f{tj)= [ f{t)dt 

n—KX) — ' JO 

J = 1 


for continuous /, e.g., the n-node composite trapezoidal rule. We will also assume that the 
quadrature rule converges uniformly over any set T of equicontinuous functions in C[0,2], 
so that if / € T , 



< e 


for n > N(e), N(e) independent of /. The n-node trapezoidal rule is an example of such a 
family, or more appropriately, the trapezoidal rule applied separately to each interval [0, 1], 


[ 1 . 2 ]. 

Nystrom’s method consists of replacing the integral in equation (4.1) with the quadrature 
rule to obtain 

- \t h (s) + J2 KMHTJt,) = g(s). (4.2) 

1 j = l 

Now let s — * * ' to obtain the 71 x 71 linear system 


- \Tn{U) + E = 9{U). (4-3) 

1 i=i 

The idea is that T n (U) « T(t x ). As shown in [2], each solution to equation (4.2) leads to 
a solution to equation (4.3) and moreover, each solution to equation (4.3) corresponds to a 
unique solution to equation (4.2) with which it agrees at the nodes fi, • ■ * ,£«. Equation (4.3) 
is the system which is solved numerically although (4.2) is the equation we will use for the 

error analysis. 

Write equations (4.1) and (4.2) as 


(-5 l + S)T = g 

(-i/ + 5„)T„ = g 


where S n is the operator in equation (4.2). Note that S n is compact since it is a finite rank 
operator on C[0,2]. Recall that S and S n depend on the parameter q. 


Theorem 4.1 T n — * T in (7(0,2] as n -+ oo. The convergence is uniform over q £ Q. 
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In proving this theorem the following result will be useful. It can be found in [2], section 3.0. 


Theorem 4.2 Let X be a Banach space, let S and (A — 5') -1 be bounded linear operators 
on X , with X ^ 0. Let T be a bounded linear operator on X with the property that either A 
is an eigenvalue ofT or (A — T) -1 exists (ifT is compact, this is satisfied). Further, assume 


ll(r-s)r|| < 

Then (A — T) -1 exists on X and 

||(A - T) _1 || < 


l(A-S , )-’lf 

l + IKA-sniim 


W-l|(A-s)-'||||(r-s)r||- 

Furthermore , if (A — S)f = g and (A — T)h = g then 


\\f-h\\<\\(\-S)-' 


H(r-5)r||||/[| + iicr-sten 
|A|-||(A- ^HIKr -S)T\\’ 


Proof of Theorem 4. 1 The proof of Theorem 4.1 is simply an application of Theorem 4.2 
with X = C[0,2], A = —1/2, S = —S and T = — S n . In the previous section it was shown 
that the map q — ► (— |/ + 5(<7)) _1 is continuous, hence || — |/ -f ^(^r)) - 1 (( can be uniformly 
bounded over Q. In order to complete the proof of Theorem 4.1 it must be shown that 


||(S’-S„)S 1 .|| - 

-4 0 

(4.4) 

ll(S--?»)9ll * 

-> 0 

(4.5) 


as n — > oo, uniformly for q £ Q. This, in conjunction with Theorem 4.2, will show that 
\\T — T n \ | — > 0 uniformly in q. 

From the argument given in section 3, for s £ [0,2] and q £ Q, the kernel K(q,s,t ) is 
uniformly equicontinuous in the t variable. The uniform convergence of ||(5' — to zero 

follows from 

{S - S n )g(s) = f (K(s, t)g(t) dt-jh K(s, t^g^j) 

Jo j = i 

< e(n) 

with e(n) — * 0 as n -4 oo, independently of q and 3. Here we have used the fact that K is 
equicontinuous in t and the assumption that the integration rule converges uniformly over 
equicontinuous sets in C7[0,2]. 
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The case for the convergence of || (S — 5 n )5 n || is similar. Let (f> be a function in (7[0,2] 
with ||</>|| = 1. Since K is uniformly equicontinuous in s and t, 

S n <f>{s + As) — S n <f>(s) — '(s + A s,tj) — K{s,tj))u)j(j>(tj) 

3 = 1 

< e(A3)p||^K| 

i=i 

where e(As) does not depend on s, t or q. For convergent integration rules the sum £" =1 \u>j\ 
is bounded in n (see [2], part I, section 4, theorem 7), so that S n (p is an equicontinuous 
function, independent of q. The rest of the argument is as in the previous paragraph but 
with g replaced by S n <f>. This shows that IKS' - S n )S n || -» 0 uniformly for q € Q and 
completes the proof of Theorem 4.1. 


4.2 Application to Inverse Problem 

Let us suppose that the temperature data for the inverse problem consists of point measure- 
ments T{ at points x; on the boundary of 1), i = 1, • • • , M . A reasonable way to attempt a 

recovery of the unknown region D is to define the functional 

M 

jm = »)(*,) - t,y 

1=1 

and seek an estimate of D(q ) as the solution to 
(IDP) minimize J(q) for q € Q 


In practice the solution to this problem will involve not the true temperature T(q), but 
the n-node Nystrom approximation T n (q). The actual minimization problem solved will be 

(IDP) n minimize J n {q) for <? G Q 


where 

M 

J"M = XKWfx,) - Tif. (4.6) 

1—1 

Based on the analysis of the convergence of T„ to T, the following theorem can be stated. 
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Theorem 4.3 Let q n be a solution to (lDP) n . Then as n -> oo, some subsequence of q n 
converges to q* £ Q. Moreover, q* is a solution to (IDP). 


Note that if the boundary data T, are consistent with some D(q*) for q* £ Q then, since 
there is a unique such q *, one can talk about the unique solution to (IDP) and a subsequence 
of q n will converge to this q*. 

Proof of Theorem 4.3 : Some subsequence of q n converges to a point q* £ Q by virtue of 
the fact that Q is compact. Let q m be any sequence in Q converging to q*. From the uniform 
convergence of T n to T over Q we can conclude that 

lim J n {q m ) - «/(<?*)• ( 4 - 7 ) 

771,71 — KX) 


For any q £ Q we have 


J n {q n ) < J n (q)- 


Taking the limit over n and using (4.7) shows that 


■W) < Aq) 


for all q £ Q, i.e., q* solves (IDP). 


5 Implementation and Examples 

5.1 Introduction 

In this section the recovery algorithm is implemented and the results of some numerical 
experiments are described. One modification is made to the assumptions of the previous 
sections: for experimental work it is more convenient to use a sample Cl which is rectangular, 
thus the boundary of Cl will not be C 2 . This make little difference to the integral equation 
formulation, for the operator (/ — S(q )), while no longer a second kind Fredholm operator, 
is in fact a small perturbation of such an operator. One can still establish the existence and 
boundedness of the inverse (/ — S(q)) see [10] for more details. The rest of the arguments 
are unchanged. The solution T on the boundary of fl is still continuous, although it will have 
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“corners” at the corresponding corners of f l. For numerical purposes it is thus beneficial to 
use a quadrature rule which allocates more nodes near these corners, rather than uniformly 
over the boundary of Q. Note that it is still assumed that dD is C 2 . 

One other modification will also be made: The heat flux g will be a point heating source, 
so that g(x) = Sp , a delta function, where P is the point at which the heat is applied. While 
this flux g does not live in the space of functions in which we have been working, as will 
be shown one can analytically remove the singularity from the problem and work with a 
smooth remainder term which fits the hypotheses made so far. In the examples that follow 
we use only the imaginary or out of phase portion of the periodic temperature response T . 
This part of the temperature is continuous, even through the singularity at the point heating 
source. The real or in phase part of T, however, has a logarithmic singularity at the point 
source. In reality one does not have a point source, but for a heat source concentrated in 
a small region the in phase response in the region of the source varies radically with the 
heat source “footprint”, which causes problems in a reconstruction algorithm. In practice 
therefore, one excludes the in phase data near the source. For convenience, we simply work 
only with the out of phase portion. 

One final remark: The linear system obtained via the n-node Nystrom’s method can be 
written 

/t(?)T„(?) = f(q) (5.1) 


where A(q) is the matrix generated by the discretization and b(q ) is the right hand side of 
the discretized integral equation; both depend on the parameter q which defines the void. 
As a result, the temperature T n also depends on q. The Levenberg-Marquardt optimization 
routine used requires the derivatives of the functional J(q) with respect to q which in turn 
requires On can obtain these derivatives by differentiating equation (5.1) with respect 
to q to obtain 


A . .8T n df dA 


(5.2) 


Both A(q) and f(q) are differentiable with respect to q and if A(q) is not singular then the 
above computation is valid, that is, ^k exists and therefore satisfies equation (5.2). Once 
T n has been obtained, ^k can be obtained from (5.2); in fact, the work done in solving the 
equation (LU decomposition) can be re-used in computing the derivative. It should be noted 
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that the temperature satisfies equation (5.2) at both the finite-dimensional (matrix) level 
and infinite dimensional (integral equation) level, i.e., also exists. 

5.2 Sample Geometry and Implementation 

The geometry used for the examples is illustrated in Figure 1. The lower left corner of the 
rectangular sample has (xi,x 2 ) coordinates (0,0). We will examine the case where the voids 
are disks of unknown center and radius, D = D(q ) with q = (a, b, r) where (a, b ) is the center 
of the disk in (xi,X 2 ) coordinates and r is the radius. The boundary of D is parameter- 
ized by xi (t) = a+r cos(27rt), x 2 (t) = b+rsm(2irt), 0 < t < 1. If the length (nq axis) of Q is L 


Point heat source 



Figure 1 : Sample geometry. 


and the height (x 2 axis) is H then the boundary of tt is parameterized from t = 1 to t = 2 


by 


(4Tt,0) 1 < f < f 


/(*) = 


(Mff(i-f)) !<*<§ 

(L-(t -§),//) l<t<l 


[ (0,H-4H(t-l)) \<t< 2. 

With appropriate bounds on the range of a, b and r, it is simple to check that this class of 
voids and parameterization satisfies the properties of section 3. The dimensions of the sample 
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for the present example will be 1.27 cm in length by 0.32 cm in height. All references to 
coordinates and sample geometry will be in centimeters. The thermal parameters correspond 
to aluminum. As mentioned above, the heat flux will be a point source applied at a points 
with (£ 1 , 22 ) coordinates (a?i, 0.32 cm) on the top of the sample. This source will have 
unit (one watt) power. Note that for the full time-dependent problem this means that the 
variation in heat flux is one watt. 

The quadrature rule used for Nystom’s method for the solution of equation (3.1) is just 
a version of the midpoint or trapezoidal rule with variably spaced nodes. Specifically, let 
f(x) = 2 p ~ l x p for 0 < x < The variable p is a parameter to adjust the spacing of the 
nodes. Allocate N nodes t{ ( N even) and weights Wi on the interval [0, 1] by 

_j/m i = i,-,f 

t{ — \ 

| 1 — f jv-»+ x i i = % + l,. . . ,N 

The corresponding weights are 

5(^1 + ^2 )> 2 = 1 

\{U+\-U- 1 ), i = 2,...,N 

1 — \{^n- 1 4- tN)i i = N. 

This rule allocates nodes and weights to (0, 1), excluding the endpoints. For p = 1 it is just 
the midpoint rule on [0, 1]. Choosing p > 1 causes the nodes to “bunch up” near 0 and 1. 
For each of the four sides of <9f l the nodes are allocated by mapping the above nodes on 
[0, 1] to the corresponding interval in t , e.g., the nodes are allocated on the bottom of fl by 
transforming the nodes on [0, 1] as t — » | + 1 with the corresponding change in the weights. 
The parameter p is chosen to be 2.0; this allocates more nodes close to the corners of fl. 
On the void boundary dD the nodes are evenly spaced with weights corresponding to the 
trapezoidal rule on [0, 1]. The typical number of nodes used is 10 to 30 on each side of fl 
and 20 to 40 on the void boundary. 

In order to deal with the singular boundary heating, first note that if T(x, y) is a funda- 
mental solution to L = (Z\ — ~) then L y V(x,y) = 0 for y ^ x, where L y means the operator 
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applied in the y variable Moreover, if P € dQ and f(y) = —2T(P,y)/a then it can be shown 
that 

ad„J(y) = 6 P -2d l/v r(P,y) on dSl, 

that is, for delta function heating T is the solution up to a nonsingular remainder term. 
The remainder term on the right hand side is actually continuous on dfi, or, more precisely, 
extends continuously through y = P. Thus one can analytically remove the singularity 
associated with the point heating and simply solve for the smoother remainder term. The 
forward heat conduction problem for the example is then 

(A - —)T = 0 in n\D (5.3) 

K 

ad.T = 2d l/y T(P,y) on d{Sl\D). 

The full solution with delta function heat flux would be given by the sum T(y ) + T(y). This 
solution has a logarithmic singularity in its real part and has smooth imaginary part. 

For reconstruction purposes, the temperature in the present examples will only be mea- 
sured along the top of the sample, so that the least squares fit functional (4.6) includes only 
this data. The heating frequency is in the range of 1 to 5 Hz. For each heating source the 
temperature response is measured at 40 equispaced points along the top of the sample. 


5.3 Strategy 

One of the necessities of an optimization approach is that one have a reasonable initial guess 
at the true void before beginning the optimization procedure, or else risk become trapped in 
a local minimum far from the “true” solution. This is particularly applicable in the present 
case, as illustrated by the following figures. The sample with void is shown in Figure 2. The 
sample is again aluminum with the same dimensions as in Figure 1. The true void D* (solid 
outline) has a radius of 0.06 cm and is centered at (£ 1 , 2 : 2 ) coordinates (0.88 cm, 0.24 cm). 
The heat source is applied directly on top of D* at 3 Hertz. The prospective void D is fixed 
to have the same radius and x-i coordinate as D*\ its x\ coordinate is allowed to vary from 
0.15 cm to 1.15 cm. The value of the functional J{q) as a function of x-i is shown in Figure 
3. The functional is of course zero when the X\ coordinates of D and D* coincide and the 
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Heating 



Figure 2: Sample geometry for least-squares functional example. 


Residual vs. void x— coordinate 



Figure 3: Functional J(q) versus D Xi coordinate. 


functional rises -steeply as one moves away from the minimum. However, if one used an 
optimization method with an initial guess which was far from the correct value {x\ < 0.5 
cm) then the optimization routine would probably not be successful, for in this region the 
functional is almost flat-actually, it slopes slightly away from the minimum. This illustrates 
the need for a reasonable initial guess at the Xj coordinate of the void D*. 

In the previous example the heat source was applied directly on top of the true void 
D*. In reality if a single heat source is applied it will not likely fall on top of or even near 
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D*. Figure 4 illustrates the same situation as Figures 2 and 3 but with the heat source 
far ( X\ coordinate 0.35 cm) from D*. Here the situation is even worse, for now the least- 
squares functional has many local minima to trap any optimization method started with 
an x i coordinate far from D*. Also, since the heat source does not “illuminate” D*, the 
functional is very flat near the minimum. In the presence of noise one would not be able to 
locate this minimum with any accuracy. 


Residual vs. void x-coordinate 



Figure 4: Functional J{q) versus D x\ coordinate. 


This leads to the following strategy for locating a void, illustrated in Figure 5: Apply 
the heat source at a number of different points along the top the sample. For each different 
heating location, take the corresponding temperature measurements along the top of the 
sample. As one passes the heating source over a void it will be detected by a change in the 
temperature response. Suppose this occurs when the heating source has an x\ coordinate 
equal to a. Then begin the optimization with initial guess x\ = a and X 2 and r anything 
reasonable. This should provide a reasonable approximation to the correct Xi coordinate of 
the void; the initial guess at X 2 and r is much less crucial. 
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Figure 5: Strategy for applying heat sources. 


5.4 Results 

This procedure is illustrated by the graphs in Figure 6. The actual void D* is located at 
(xi,X 2 ) coordinates (0.77 cm, 0.24 cm) with radius 0.06 cm. A total of 9 point heating 
sources at 3 Hertz were applied on the sample top surface, equispaced, spanning one half 
of the sample length, with source 5 centered on the top surface. The graphs show the out 
of phase or imaginary portion of the temperature response for each source location. The 
temperature response changes most rapidly and peaks between heating locations 3 and 4. 
Contrast this to the temperature response when no void is present, shown in Figure 7; the 
response does not change, since nothing under the heating source changes as the source 
moves. A Levenberg-Marquardt algorithm implemented along the lines in [8] was used to 
recover an estimate of the void D *, using only the peak temperature response data, heating 
sources 3 and 4. The initial guess at the coordinate was halfway between these sources 
at 0.75 cm. The x 2 and radius initial guesses were 0.16 cm and 0.1 cm, respectively. The 
optimization code converges to the correct void in 1 1 iterations, reducing the residual from 
0.154 to less than 0.002. However, no noise was present. 

As a more realistic example, we take the same geometry and void as the previous ex- 
ample but with 20 percent zero-mean gaussian noise (measure as percent of signal RMS 
value) added to the temperature response. The responses are shown in Figure 8. The largest 
response is for position 3. The data for heating positions 2, 3 and 4 were then used in the 
optimization with initial guess x\ = 0.79 cm (the X\ coordinate for source 3), x 2 — 0.16 cm 
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and 7 * — 0.1 cm. The initial residual is 0.406. The optimization routine reduces this to 0.381 
in 12 iterations. The final estimate for the void is x\ = 0.755 cm, x 2 = 0.22 cm, r = 0.067 
cm. The true void is shown as the solid outline in Figure 9 and the and recovered estimate 
as the dotted outline. 


Pos. 1 


Pos. 2 


Pos. 3 



Pos. 4 



Pos. 7 




Pos. 5 



Pos. 8 




Pos. 6 



Top of sample, cm 
Pos. 9 



Top of sample, cm 

Figure 6: Out of phase temperature for varying heat source locations. 
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Figure 7: Out of phase temperature for varying heat source locations, no void. 
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Figure 8: Out of phase temperature for varying heat source locations, 20 percent noise. 
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6 Conclusion 


The problem of recovery a void in a material sample based on the sample’s surface tem- 
perature response to external heating has been considered. Uniqueness and continuous 
dependence results have been shown, and an optimization based algorithm for recovering 
an estimate of the void has been demonstrated. Much of this work rests on reformulating 
the heat conduction problem as a boundary integral equation, which provides a means of 
rapidly solving the heat conduction equations. The algorithm was run for the simple case 
of a circular void and computationally generated data. This algorithm exhibits some of the 
problems that optimization approaches are heir to, specifically, the likelihood that a poor 
initial guess will not converge to a global minimum, and some strategies for overcoming this 
have been described. 

Collection of actual data from a experimental setup at NASA Langley Research Center 
has already been performed. Preliminary analysis of this data shows good agreement with 
the heat conduction model and the ability to actually recover subsurface voids. The analysis 
of this experimental data will be reported elsewhere. Of interest for future research is a study 
of the sensitivity of this thermal technique, e.g., the size of voids which can be detected at 
various depths. Techniques for voids of different shapes (especially cracks or disbonds) should 
also be examined. The heat conduction model could also be improved; zero boundary flux 
away from the source becomes unrealistic at low frequencies. We would also like to pursue a 
full three-dimensional heat conduction model leading to a two-dimensional boundary integral 
formulation. Such a formulation would require a finite or boundary element technique for 
the numerical solution, rather than Nystrom’s method. 
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